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This is the final paper in a series that introduces geodesic molecular dynamics at constant po- 
tential energy. This dynamics is entitled NVU dynamics in analogy to standard energy-conserving 
Newtonian NVE dynamics. In the first two papers [Ingebrigtsen et at, J. Chem. Phys. 135, 104101 
(2011); ibid, 104102 (2011)], a numerical algorithm for simulating geodesic motion of atomic systems 
1 was developed and tested against standard algorithms. The conclusion was that the NVU algorithm 

has the same desirable properties as the Verlet algorithm for Newtonian NVE dynamics, i.e., it is 
. time-reversible and symplectic. Additionally, it was concluded that NVU dynamics becomes equiv- 

' alent to NVE dynamics in the thermodynamic limit. In this paper, the NVU algorithm for atomic 

systems is extended to be able to simulate geodesic motion of molecules at constant potential energy. 
We derive an algorithm for simulating rigid bonds and test this algorithm on three different systems: 
an asymmetric dumbbell model, Lewis- Wahnstrom OTP, and rigid SPC/E water. The rigid bonds 
introduce additional constraints beyond that of constant potential energy for atomic systems. The 
CNj , rigid-bond NVU algorithm conserves potential energy, bond lengths, and step length for indefinitely 

long runs. The quantities probed in simulations give results identical to those of Nose-Hoover NVT 
dynamics. Since Nose-Hoover NVT dynamics is known to give results equivalent to those of NVE 
dynamics, the latter results show that NVU dynamics becomes equivalent to NVE dynamics in the 
thermodynamic limit also for molecular systems. 
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I. INTRODUCTION 



Iii two recent papers*^ (henceforth: Papers I and II) molecular dynamics at constant potential energy was intro- 
O ■ duced, tested, and compared to well-known molecular dynamics algorithms. This new molecular dynamics is entitled 
NVU dynamics in analogy to standard energy-conserving Newtonian NVE dynamics. The conclusion was that NVU 
dynamics is a fully valid molecular dynamics, which for sufficiently large systems can be used interchangeably with 
7— I ■ NVE dynamics for calculating most quantities of interest. NVU dynamics is not faster than standard NVE or NVT 
dynamics, but introduces a new way of thinking about molecular dynamics. Molecular dynamics at constant potential 
,— ' ■ energy was previously considered by Cotterill and co-workers^—, by Scala et alX, and most recently by Stratt and 
co-workers^—, who actually allowed also lower potential energy values. Our motivation for studying NVU dynamics 
derived from recent work on strongly correlating liquids and their isomorphsi^— (see the Introduction of Paper I). 
NVU dynamics is defined by geodesic motion on the constant-potential-energy hypersurface £1 defined by 

OV 

o: 

CN ■ O _ fD r- U3N 



!l = {Re i? JiV | C/(R) = U }. (1) 

Here R = {rW, r^'} in which rW is the position vector of the k'te particle (we follow here the notation of the 
Appendix of Paper II) , and U is the potential-energy function of an A^-particle classical system. A geodesic on fl is 
a curve that satisfies the condition of stationary length for fixed endpoints Ha and Rg, i.e., 



6 dl 
Jr a 



-0, (2) 
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where dl is the line element of the metric. The shortest path between any two points is a geodesic. On a sphere 
geodesies are great circles, the "straightest lines" of the surface. Traversing a geodesic at constant velocity thus 
corresponds to a generalization of Newton's first law to a curved space (the surface itself). 

In Paper I the NVU algorithm was developed via a discretization of Eq. ([2]), subsequently carrying out the variation. 
This technique, which is known as variational integration^—, resulted in a "basic" NVU algorithm that is similar to 
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the well-known Verlet algorithm Rj+i = 2R; — Rj_i + (At) 2 Fi/m for Newtonian (NVE) dynamics (m is the particle 
mass which is assumed identical in this section, and Fj = -Vr ; (7 is the 3-ZV-dimensional force vector); the index i 
refers to step i of the integration sequence. In the Verlet algorithm At is a fixed time step length. In comparison, the 
basic NVU algorithm is given by (Paper I) 



R 



2Ri — Ri_i + 



-2F, • (Ri — Ri-i) 



(3) 



If the number of particles N increases, the relative variation of the term — 2Fi • (Ri — Rj_i)/F? decreases, and 
this is why equivalence with Newtonian NVE dynamics is established in the thermodynamic limit. This equivalence 
should be understood in the sense that the relative deviations between, for instance, NVE and NVU auto-correlation 
functions go to zero as N —> oo. 

Paper I additionally developed a "stabilized" version of the basic NVU algorithm to prevent accumulation of 
numerical errors. This version of the algorithm is given by (defining the position changes A i+1 / 2 = Rj+i — R;) 



where Iq is the step length and 
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lAi-i/all' 
Rj+i = R t + A i+1 / 2 , 



(4) 
(5) 
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(-2F, • Ai_ 1/2 + Ui-x - U ) 



(6) 



All simulations in Papers I and II were performed with the stabilized algorithm. The basic algorithm was used, 
however, for theoretical considerations. 

In this article we extend the stabilized NVU algorithm to deal with simulations of molecular systems. Molecular 
systems are simulated by introducing rigid and/or flexible bonds between the atoms in the modelling. Flexible bonds 
introduce merely an additional contribution to U, for instance, harmonic spring potentials. The basic NVU algorithm 
conserves the total potential energy and can readily simulate flexible bonds. The focus in this paper is thus on 
implementing rigid bonds in the framework of NVU dynamics. 

Sect ion ITT1 considers NVU dynamics with rigid bonds. Introducing rigid bonds in the simulations leads to Lagrangian 
multipliers in addition to those introduced in order to keep the potential energy constant (Paper I). Section [TT] is fairly 
technical and easiest to read after reading Paper I. Section UTT] gives simulation and model details. Section HVl tests 
the rigid-bond NVU algorithm, and Sec. Ivl investigates the NVU sampling properties by comparing the NVU results 
to Nose-Hoover NVT result o 24 ' 25 on three different systems: the asymmetric dumbbell model 2 ^, Lewis- Wahnstrom 
OTP 2 —, and rigid SPC/E water—. Nose- Hoover NVT dynamics is known to give results equivalent to NVE dynamics 
in the thermodynamic limit—, and we refer to these dynamics interchangeably in the forthcoming sections. Finally, 
Sec. rVTl concludes. 

II. RIGID-BOND NVU ALGORITHM 

The rigid bond s 30 ' 31 introduce constraints among the particle coordinates of the system. Each constraint a = 1, G 
is of the form 



<r Q (R) = (r< fc «> - r^)) 2 = (r Q ) 2 = C 2 a , 



(7) 



it expresses that the distance between particles k a and l a is a constant, C a . In Papers I and II the integral of Eq. 
([2]) was merely restricted to the constant-potential-energy hypersurface fl. Each rigid bond constraint introduces a 
function a a to be kept constant, and thus the integral of Eq. ^fy is now further restricted to the sub-manifold u) of VL 
where the bond constraints are satisfied, 



{Refi & a (R) 
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1,...,G}. 



(8) 
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If the bond constraints are independent, as assumed throughout the paper, a; is a (3N - G - l)-dimensional compact 
Riemannian manifold. The variational principle defining NVU dynamics with rigid bonds is given by 



5 dl = 
Jr a 



(9) 



Most of Papers I and II dealt with the case of identical particle masses, but we wish here to develop a completely 
general molecular NVU algorithm. The line element dl is defined by 



d^5> fc (dr«) 2 , (10) 

k 

where fhu = mk/(m) is the "reduced" mass of particle fc. Equation (|10[) is not the standard Euclidean line element, 
but a mass- weighted line element that goes back to Hert a 32 i 33 . We shall refer to this metric as the "Hertzian" metric. 
The point of this particular metric is that it ensures equivalence between NVU and NVE dynamics for systems of 
atoms and molecules of varying mass. In appendix [X] we derive the variable-mass atomic NVU algorithm applying 
the Hertzian metric (correcting also a typo of the Appendix of Paper II). 
Applying the variational integration technique to Eq. (|9]) gives 



6 E wE ™* ( r * (fe) - r <-i) 2 - E w R ') + E A«i*aCRi) = o . (ii) 

\ i \ k i i,a J 

In Eq. (|11[) the path is divided into a number of discrete points and one Lagrangian multiplier A a i is introduced for 
each constraint a at every point i. Following standard notation for constraint molecular dynamic a 30 ' 31 , the Lagrangian 
multipliers of the bond constraints are chosen with a positive sign. As in Papers I and II we now make the Ansatz of 
constant step length l , i.e., 

Em^rf-r^) 2 ^. (12) 

k 

Carrying out the variation of Eq. (fTTj) using Eq. (|12|) leads to (compare the derivation in Paper I) 



rffi = 2r| fe) - r& + J^ff > + A. V ^ £ A m a a , (13) 

lilfc ink 1 

a 

where = — V mU is the force on particle k at step i. This equation constitutes the NVU algorithm with rigid 

bonds. It has a close resemblance to the Lagrangian equations of motion with holonomic constraints 3 ^, i.e., rigid- 
bond NVE dynamics^. Equation ([T3j) contains G + 1 Lagrangian multipliers for each integration step, which must 
be determined to complete the algorithm. 



A. Determining the NVU Lagrangian multipliers 

This section shows how to calculate the Lagrangian multipliers. Since the algorithm is to be implemented on a 
computer (with finite-precision), we shall proceed directly to a "stabilized" algorithm conserving for indefinitely long 
runs potential energy, bond lengths, and step length (in 3A-dimensions). The resulting algorithm reduces to the 
stabilized atomic NVU algorithm of Eqs. (H])-© in the case of no bonds constraints. 

Some notation used in the following derivation is now introduced (the nomenclature of text is summarized in Table 
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Symbol 


Definition 






m fe 


The a'te bond constraint between particles fc a and Z a with a = 1, G. (a a = (r°) 2 = C^). 
The mass of particle divided by the average mass of the system, (rhk = m k /(m)). 




3-dimensional vectors 


°i+l/2 
j(fe) 

* W 

r a 

1 % 

SOL 

°i-l/2 

r 

% 

Si 


Position of particle k at step i. 

Displacement of the position of particle k between step i and i + , 2 = r i+i ~ r i*^)- 

Force on particle k at step i. (f^' = —V (k)U). 

Constraint force on particle k at step i. (g^ = V (k) A a ;<7 Q ). 

Displacement of the positions of particles k a and l a at step i. (rf = r ( fc °) _ r£ ). 

Displacement of the velocities of particles k a and l a at step i — 1/2. (<5f_i/ 2 = ^j*™^ — ^i-i/2) 1 

Sum of displacements of positions and velocities of particles k a and Z a at, respectively, step i and 

z-1/2. (sf = rf + «f_ x/a ). 

Displacement of the forces on particles fc a and Z a at step i divided by their reduced particle mass. 

Displacement of the constraint forces on particles k a and l a at step i divided by their reduced 
particle mass, (g" = g[ ka) /rh ka - gV /rhi a ). 




3-/V-dimensional vectors 


Ri 

F, 
F, 
G, 


Position of all particles at step i. (Hi = {r' 1 ', ...,r^ }). 

Displacement of the positions between step i and i + 1. (Aj+1/2 = R;+i — Ri). 
Force on all particles at step i. (Fj = — Vp^f/). 

Force on all particles at step i divided by the reduced particle mass. (Fj = {f} 1 '/mi, fj /mjv}). 
Constraint force on all particles at step i divided by the reduced particle mass. (Gj 

, (1) / - (JV) / - IX 

{g, 7mi,...,gJ Vmjv}). 



TABLE I. Definitions and nomenclature of the text. 



Defining 5| +1 / 2 = r^ + \— r s - and g ■ = V r (*o J2 a ^ aiCJa the "Leap-frog"— version of the rigid-bond 7VT^J7 algorithm 
Eq. reads 



C\/ 2 = e\/ 2 + -A#> + A g « (14 ) 

In analogy to rigid-bond NVE dynamics we call the "constraint force" on particle k at step i. Introducing 

the notation F l = {fj^/mi, /rh N } and G, = {gf 3 /"id, g,^ /m N }, the algorithm in the full 3N- 

dimensional coordinate space reads 



A l+1/2 = Ai_ 1/2 + ioA-jFj + ZoGi, (16) 
R m =R, + Ai+i/a,. (17) 

The Lagrangian multipliers are calculated by combining a result derived in Paper I with the method applied in the 
SHAKE algorithm 30 for rigid bonds in NVE dynamic s 30 ' 35 ' 36 . The SHAKE algorithm calculates the Lagrangian 
multipliers from the equations (rf +1 ) 2 = C\. In doing so, the target value of the constraints C a appears explicitly 
in the algorithm, making the bond lengths insensitive to numerical error. The expression for r" +1 is supplied by the 
integration algorithm containing herein the Lagrangian multipliers. In our case, this gives G equations with G + 1 
unknowns. The missing equation is supplied by an expression derived in Paper I, namely that C/^+i = t/$_i — Fj • 
(Ri+i — Rj-i) to third order in the step length. In the discrete sequence of points f/j+i is set equal to Uq (the constant 
defining f2), making the constraint of constant potential energy also insensitive to numerical error. We thus have the 
following G + 1 equations for calculating the Lagrangian multipliers 
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Ui-i - F, • (R i+1 - Rj_i) - f/ = 0, 

(if +1 ) 2 -C 2 =0, (a = l,...,G). 

ByEqs. (02} and (TU); Rj +i -Rj i = A i+1/2 + Aj_ 1/2 = 2Aj_ 1/2 + / A i F 1 + / Gj. Defining 



(18) 
(19) 



1/2 



t(*a) 



ff = ^°V™fc„ - H^/micJ and 9? = Si/r^kc, - gf a) M/ Q > since b Y Ec l s - and G3>; r f+i = r i+i - r< i+i 



i-1/2 "i-1/2' 

(fee) „(ia 



S ?+i/2 - *&/a = r ? + 5 ?-V2 + ^AiC + fog?, it follows that 



i-l 



Fi • 2A l _ 1/2 + foAjFj + Z Gi| - 17 = 0, 

2 



0, (a = l,..,G). 



(20) 
(21) 



The above coupled quadratic equations for the Lagrangian multipliers are now solved following the produce of the 
MILC-SHAKE algorithm 37 , which starts by neglecting the second order terms in the Lagrangian multipliers and 
solving the resulting linear equations. Afterwards, the second order terms are taken into account in an iterative 
manner - the details of which are described below. 

For each integration step i, the linearized equations are given by 



A,;Ai 



(22) 



where Aj is a (G + 1) x (G + 1) matrix, Aj = {Aj, Ai,-, Aci}, and bj a G + 1 column vector. We start by calculating 
explicitly the first few elements of the matrix Aj. An consists merely of the factor infront of Aj in Eq. (|20l) . 
i.e., An = —lo^i ■ Fi. The second element A12 appears after expansion of the dot product Fj • Gj. Noting that 
V r ( fcQ ,a Q = 2if , we have F, G, = ff> • g^/mi + - + tf • sf ] /rh N = 2A H (f- • rj) + ... + 2A Gl (ff • rf ). The last 
equation follows as the Lagrangian multipliers appear in pairs, differing only by the sign from V (* a )Oa and the term 

f\ kc "^ /fhk a ■ We thus find A12 = —210% ■ r ji ^13 = — 2Zoff ■ rf, etc. In the second row of Aj, the short-hand notation 
sf = rf +5j*_]/2 is introduced, making A21 — 2lo( s } -%), i- e -, the factor infront of Aj after squaring of the parentheses. 
The next element A22 appears after expanding sj • — sj ■ £^ A ( 3j(^i— V r (fci)<r | a — :pj~V r Oi)crg). In this sum, we 
identify the factor in front of Aij, giving A22 — 2losj ■ (-J— V 00 ci — V and similarly for the remaining 

^1 r i 'l r i 

elements of the second row. 

Altogether, the elements of Aj are thus given by 



/-F< • F./2 



2Zo 



-f G • r G 



s} • E 



S, 1 • (-^-V ( fcl )CTi - J-V (lOITl) 



S,- • (^^V (k^CTn (Zi)(Jg) 



The column vector consists of all zeroth-order terms in Eqs. ([20]) and ([21 



(23) 



/[/o-f/i-i+2Fj- Aj_ 1/2 \ 
C\ - (si) 2 



V 



r 2 



(sf) 2 



/ 



(24) 



Turning now to the iteration procedure, the second-order terms in the Lagrangian multipliers (Eq. ([21]) ) are taken 
into account by iterating the right-hand side of Eq. (f2"2"j) via the scheme (a — 1, G) 



+1 _ 



Cl ((r? +1 ) 2 ) J 



(25) 
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The superscript j refers here to iteration j, and ((r" +1 ) 2 )' J are the positions associated with iteration j. The element 
60 is not updated as it derives from the constraint of constant potential energy. For each iteration j the term 
Cq — (( r f+i) 2 )' 7 is expected to become smaller as the bonds are satisfied better and better, and convergence was 
achieved within a few iterations^!. 

For each integration step i, the algorithm for determining the NVU Lagrangian multipliers thus proceeds as follows 

1. The Lagrangian multipliers of iteration j, (Aj) J , are calculated from Eq. (|22p. 

2. ((rf +1 ) 2 ) 3 is calculated via Eqs. flU) and JH]) using (Aj) J . 

3. hi is updated via Eq. (25]) from ((rf +1 ) 2 ) J . 

4. The above steps are repeated (starting iteration j + 1) until convergence is established (we used a preset number 
of iterations, typically 3-5). 

How is constant step length Iq ensured numerically? Generalizing the approach of Paper I we introduce a normalizing 
factor such that 

(fc) 

1+1/2 "^nn' ( } 

r (fc) _ r (*) + *(k) (27) 

r i+l r i + °i+l/2> V Z 'J 

where 

The normalizing factor is close to unity-^ and ensures trivially J2k ^ki^^^Y = ^Q) i- e -> that the step length is 
conserved. The algorithm is now absolutely stable, conserving potential energy, bond lengths, and step length for 
indefinitely long runs. The stability of the NVU algorithm is tested numerically in Sec. IIV1 



B. Alternative determination of the NVU Lagrangian multipliers 

The previous section followed the traditional way of calculating the Lagrangian multipliers. The NVU Lagrangian 
multipliers may also be calculated by Taylor expanding the constraints a a in analogy to the method sketched above for 
the potential energy. In this way, the constraints of constant potential energy and constant bond lengths are treated 
on equal footing. The set of equations to be solved is the following (recall that R,+i R; 1 = 2A i _ 1 / 2 + ^o^iFi-|-/oGi) 

^-1 - F t ■ (Ri+i - Ri_i) - U = 0, (29) 
ffa( l -i)+VR ! a m -(R l+1 -R l _i)-C 2 =0, (a=l,...,G). (30) 

The determination of the Lagrangian multipliers is linear and thus no iterations are needed. The bond constraints 
a a are obeyed to the same order 0(Iq) as the constraint of constant potential energy. The sampling properties 
of this novel, alternative determination method is tested briefly in Sec. |Vl It appears to be a promising new way 
of determining the Lagrangian multipliers in connection with rigid bonds, which might also be useful for standard 
bond-constraint NVE or NVT simulations. 



III. SIMULATION DETAILS AND MODEL SYSTEMS 



We investigated three systems: The asymmetric dumbbell model, the Lewis- Wahnstrom OTP model, and rigid 
SPC/E water. For all simulated pair potentials the shifted-force truncation scheme was applied at a cut-off radius r c . 
If the pair potential is v(r) and the pair force is f(r) = —v'(r), the shifted force is given 
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Mr)= ' (r) - /(Pc) l[ r<rC ' (3D 
10 it r > r c . 

This corresponds to using the following pair potential below r c : v$p{r) = v(r) — v'(r c )(r — r c ) — v{r c ). All simu- 
lations were performed with the NVT and NVU algorithms. Recall that NVE and NVT dynamics give equivalent 
results 2 ^; for this reason no simulations are presented for NVE dynamics. The RUMD code^S was used for molecular 
dynamics simulations (an optimized open-source GPU code). The NVT ensemble is generated via the Nose- Hoover 
algorithm— i 25 ' 40 , and the bonds held fixed using the time-reversible constraint algorithm of Refs.— and — . The NVU 
algorithm is described in Sec. [TTJ The starting files for NVU dynamics were taken from an equilibrated NVT simu- 
lation. The positions and velocities of the NVT configuration do not correspond perfectly to motion on uj, since the 
potential energy and step length are not those of Uq and lo, respectively. As all three constraints are to be satisfied 
simultaneously, this results in numerical problems when starting the simulation from the particular NVT configura- 
tion. A more gentle procedure is thus applied, where the atomic NVU algorithm is used for a couple of integration 
steps to ensure the values of Uq and Iq. Afterwards, the rigid-bond NVU algorithm is used. 



A. NVU iteration procedure 

The quadratic equations (Eq. (|25j) ) were iterated with a fixed number of iterations (between 3 and 5). The linear 
systems were solved utilizing CUSP 41 , a library for solving systems of linear equations on the GPU. More specifically, 
the stabilized biconjugate gradient algorithm with a Jacobi preconditioner 4 — was used with the initial value Xi = 0. 
The relative tolerance r for the asymmetric dumbbell and Lewis- Wahnstrom OTP models was chosen as r = 10~ 7 
and for rigid SPC/E water as r = 3 • 10~ 7 . A larger tolerance was chosen for rigid SPC/E water due to convergence 
issues in connection with shifted-force Coulomb interactions (see below). 

The maximum number of allowed iterations was 50. A restart scheme was applied when the solver did not converge 
within the chosen tolerance. In this case the solver (and quadratic iteration) was restarted from the partially estimated 
"solution" adding 2 • 10 -7 to the tolerance. It should be noted that the stabilized biconjugate gradient algorithm may 
get "trapped", resulting in a break-down of the CUSP linear solver. If this happens, it is detected by our program, 
and the solver and quadratic iteration are restarted, with a smaller number (10) of maximum allowed iterations for 
the solver. 



B. The asymmetric dumbbell 

The asymmetric dumbbell model^ 6 . consists of a large (A) and a small (B) Lennard- Jones (LJ) particle, rigidly 
bonded with bond distance of tab = 0.29/0.4963 (here and henceforth units are given in LJ units referring to the 
A particle such that <jaa — U £AA = 1, and ttia = 1)- The asymmetric dumbbell model has ctbb = 0.3910/0.4963, 
(bb = 0.66944/5.726, and tub = 15.035/77.106. The AB interaction between different molecules is determined by 
the Lorentz-Berthelot mixing rule 34 , n — 500 molecules (here and henceforth n denotes the number of molecules and 
N the number of atoms) were used in the simulations with a pair-potential cut-off of r c — 2.5. The step length l was 
fixed in the range 0.125-0.138 depending on the state point. 

Simulations were also performed where the rigid bonds were replaced with stiff harmonic springs. The spring 
constant was k = 3000, while all other model parameters remained unchanged. 



C. Lewis- Wahnstrom OTP 



The Lewis- Wahnstrom OTP model 27 consists of three identical LJ particles rigidly bonded in an isosceles triangle 
with sides of taa = 1 and top angle of 75°. All parameters (including the masses) are unity for the OTP model, n 
= 320 molecules were simulated and a pair-potential cut-off of r c = 2.5 was used. The step length was 0.100. 



D. SPC/E Water 

The SPC/E water model 2 ^ is an isosceles triangle with sides toh — 1/3.166 and top angle 109.47°. The OO 
intermolecular interactions are given by the LJ pair potential (eoo = L &oo — L and mo — 15.9994/1.00794). The 
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three particles are charged with qo = —22.0 and qn = \qo\/2- n — 2000 molecules were simulated and a pair-potential 
cut-off of r c = 6.28 for both LJ and Coulomb interactions was applied 43,44 . The step length was fixed in the range 
0.06-0.07 depending on the state point. For this system the numerical stability is surprisingly sensitive to the cut-off 
used in the Coulomb interactions, but a larger shifted-force cut-off improves this behavior—. 

IV. TESTING THE STABILITY OF THE RIGID-BOND NVU ALGORITHM 

This section tests the conservation properties of the rigid-bond NVU algorithm. Table [TT| shows the potential 
energy, the deviation of bond lengths, and step length as functions of integration step number for Lewis- Wahnstrom 
OTP at p = 0.329 and T = 0.700. It is clear that these quantities are conserved by the algorithm and that no drift 
occurs. The step length is conserved to the highest accuracy since it is not prone to numerical error in determining 
the Lagrangian multipliers. 



Integration steps 


U 1 N 


(1/G£ a (r a -C a ) 2 ) 1/2 




10 1 


-4.42550 


2.81207 10~ 7 


0.0999999 


10 2 


-4.42552 


3.03535 10~ 7 


0.1000000 


10 3 


-4.42552 


2.81128 10" 7 


0.1000000 


10 4 


-4.42552 


2.95078 10" 7 


0.1000000 


10 5 


-4.42550 


3.08793 -10" 7 


0.1000000 


10 6 


-4.42551 


2.90477 10" 7 


0.1000000 



TABLE II. Potential energy, deviation of bond lengths and step length as functions of integration step number in the NVU 
algorithm for Lewis- Wahnstrom OTP (p = 0.329, T = 0.700). Single-precision floating-point arithmetic was used for the 
simulations. 

Figure [T] shows the distribution of the term lo\i(m) in Eq. (|13[) (recall rhk = rrik/(m}). In NVU dynamics there 
is, as such, no notation of time; a geodesic on the manifold can be traversed with any velocity. Comparing the NVU 
algorithm of Eq. (fL3"|) to the rigid-bond Verlet algorithm^ r-^ = 2r^ — + ((^) 2 / m k)[^ +g[ k ^], we can define 
the term l \ i {m) as a varying "time step" length of the NVU algorithm (see also Paper II), i.e., 

{At i<NVU ) 2 =l Xi(m). (32) 

The integration steps of the NVU algorithm are thus henceforth referred to as "time steps". The average of Eq. 
(|3"2")l is used in Sec. IVl when comparing to NVT dynamics. As was the case for the atomic NVU algorithm (Paper 
I), loXi(m) is Gaussian distributed for large systems and its relative variation decreases as the number of particles 
increases. It thus becomes a better and better approximation to treat this term as constant, implying equivalent 
sampling properties of NVU and NVE dynamics also when rigid bonds are included in the simulations. 



lxio 6 rr 




lA.<m> 

i 

FIG. 1. The probability density of the "time step" length (Ati,NVu) 2 = ioAi(m) of the rigid-bond NVU algorithm for Lewis- 
Wahnstrom OTP at p = 0.329 and T = 0.700. n = 320 molecules were used in the simulations. 
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V. SAMPLING PROPERTIES OF THE RIGID-BOND NVU ALGORITHM 

The NVU algorithm is now compared to NVT dynamics for the three different models. First, we consider the 
asymmetric dumbbell model^ 6 -, both rigid and flexible. Afterwards, the Lewis- Wahnstrom OTP models, and finally 
rigid SPC/E water-22. 



A. The asymmetric dumbbell model 

In Figs. Ufa) and (b) are shown, respectively, the molecular center-of-mass (CM) radial distribution functions 
and the CM incoherent intermediate scattering functions for the rigid asymmetric dumbbell model 26 for different 
temperatures at p = 0.932. The black circles and curves give NVT simulation results while the red crosses give the 
NVU simulation results. The two radial distribution functions in Fig. [^a) agree very well, and this is also the case 
for the dynamics in Fig. EJb). 




FIG. 2. Comparison of structure and dynamics in NVU and NVT simulations of the rigid asymmetric dumbbell model. The 
black circles and curves give NVT, the red crosses NVU simulation results, (a) The molecular CM radial distribution functions 
at p — 0.932 and T = 0.500. (b) The molecular CM incoherent intermediate scattering functions at p = 0.932 and T = 0.500, 
0.600, 0.700, 0.800, 0.900. 

For reference, we also simulated (Fig. [3]) the corresponding quantities for the flexible-bond asymmetric dumbbell 
model at the state points of Fig. [2] Again, there is a very good agreement between NVU and NVT dynamics. 




FIG. 3. Comparison of structure and dynamics in NVU and NVT simulations of the flexible-bond asymmetric dumbbell 
model. The black circles and curves give NVT, the red crosses NVU simulation results. The same state points as in Fig. [2]were 
simulated, (a) The molecular CM radial distribution functions at p — 0.932 and T = 0.500. (b) The molecular CM incoherent 
intermediate scattering functions at p = 0.932 and T = 0.500, 0.600, 0.700, 0.800, 0.900. 
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B. Lewis- Wahnstrom OTP 



We show in Figs. Hla) and (b), respectively, the molecular CM radial distribution functions and CM incoherent 
intermediate scattering functions for the Lewis- Wahnstrom OTP model 27 . The same symbols and meanings as in the 
preceding section are used. Again, the NVU and NVT simulations agree very well for both structure and dynamics. 




FIG. 4. Comparison of center-of-mass structure and dynamics in NVU and NVT simulations of the Lewis- Wahnstrom OTP 
model. The black circles and curves give NVT, the red crosses NVU simulation results, (a) The molecular CM radial distribution 
functions at p — 0.329 and T = 0.700. (b) The molecular CM incoherent intermediate scattering functions at p — 0.329 and T 
= 0.700, 0.800, 0.900, 1.000. 



For comparison, we also show in Fig. [5] the corresponding particle quantities for the OTP model. 




FIG. 5. Comparison of particle structure and dynamics in NVU and NVT simulations for the Lewis- Wahnstrom OTP model. 
The black circles and curves give NVT, the red crosses NVU simulation results, (a) The particle radial distribution function 
at p — 0.329 and T = 0.700. (b) The particle incoherent intermediate scattering functions at p = 0.329 and T = 0.700, 0.800, 
0.900, 1.000. 



C. SPC/E Water 



Finally, we consider in Fig. [6] the same quantities as above for the rigid SPC/E water model 28 . Again, full 
equivalence between NVU and NVT dynamics is found. 
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FIG. 6. Comparison of structure and dynamics in NVU and NVT simulations of rigid SPC/E water. The black circles and 
curves give NVT, the red crosses NVU simulation results, (a) The molecular CM radial distribution functions at p — 1.000 
and T = 3.800. (b) The molecular CM incoherent intermediate scattering functions at p = 1.000 and T = 3.800, 4.200, 5.000. 



The linear algorithm for determining the Lagrangian multipliers presented in Sec. Ill Bl (Eqs. ([25)1 and ([50)) ) is 
tested in Fig. [7] by probing the molecular CM radial distribution functions. NVU and NVT dynamics also here give 
identical results. 




FIG. 7. Comparison of structure in NVU and NVT simulations of rigid SPC/E water at p = 1.000 and T = 3.800 applying 
the linear method to determine the Lagrangian multipliers (Eqs. (|29p and (|30p h The bond lengths are here conserved to order 
10~ 6 in the standard deviation of the bonds (using single-precision). 

We conclude from the presented results that for sufficiently large molecular systems with flexible and/or rigid bonds, 
NVU dynamics is equivalent to Nose- Hoover NVT dynamics (and, by implication, Newtonian NVE dynamics). 



VI. SUMMARY 



NVU dynamics is molecular dynamics at constant potential energy realized by tracing out a geodesic on the 
constant-potential-energy hypersurface 51 (Eq. (JTJ). In Papers I and II-- 1 -, a "basic" and a "stabilized" atomic NVU 
algorithm for simulating geodesies on f2 was developed. The basic NVU algorithm has excellent stability; it is time- 
reversible and symplectic, and the stabilized algorithm was developed only to prevent accumulation of numerical error. 
It was found that atomic NVU dynamics becomes equivalent to atomic NVE dynamics in the thermodynamic limit. 

In this paper the stabilized NVU algorithm has been extended to simulate molecules at constant potential energy. 
Molecules are simulated by introducing rigid and/or flexible bonds in the models. The atomic NVU algorithm keeps 
the potential energy constant and can thus right away simulate flexible bonds. The focus here was on incorporating 
rigid bonds in the framework of NVU dynamics, which leads to the introduction of additional Lagrangian multipliers 
beyond those of the constraint of constant potential energy. This is completely analogous to the approach for simu- 
lating rigid bonds in standard Newtonian NVE dynamic a 30 ' 35 i 36 . In the NVU algorithm, a set of coupled quadratic 
equations was constructed for calculating the Lagrangian multipliers and solved in an iterative manner as a linear 
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system, a procedure developed for rigid-bond NVE dynamics in the MILC-SHAKE algorithm 37 . In addition, a set 
of linear equations was presented for calculating the Lagrangian multipliers, and appears to be a promising new way 
of simulating rigid bonds. 

The rigid-bond NVU algorithm reduces to the atomic NVU algorithm when there are no rigid bonds. The algo- 
rithm was tested on three different model systems: the asymmetric dumbbell model, Lewis- Wahnstrom OTP, and 
rigid SPC/E water. The probed quantities in the simulation gave identical results to those of Nose- Hoover NVT 
dynamics. We conclude that also for molecular systems do NVU dynamics become equivalent to NVE dynamics in 
the thermodynamic limit (since NVE and NVT dynamics are known to give equivalent results 2 ^). 
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Appendix A: Derivation of the atomic NVU algorithm for the Hertzian metric 

According to Newtonian dynamics, heavy particles move slower than light particles in thermal equilibrium. The 
standard Euclidean metric does not involve the particle masses, and thus applying this metric to geodesic motion for 
systems of varying masses will not produce dynamics equivalent to Newtonian dynamics in a thermal system. The 
mass-weighted metric of Herta22, however, ensures that NVU dynamics becomes equivalent to NVE dynamics in the 
thermodynamic limit, as is clear from the derivation below. This metric is given by (where rhk — mk/(m)) 

dl 2 =£m fc (drW) 2 ! . (Al) 

fc 

We here derive the discrete NVU algorithm applying the Hertzian metric (this appendix also corrects a typo in Eq. 
(A5) of Paper II). The discretized variational condition for geodesic motion on f2 is 

6 (E ^£ A * (r^-ri^) 2 - E W R *)) = " ( A2 ) 
Assuming a constant step length Zq, i.e., 

E^(^-^i) 2 -^ (A3) 

fc 

it follows by differentiation with respect to from Eq. (|A2[) that 

m*(rf - r W) + m fc (rf > - r«) + Z A 4 ff } = 0. (A4) 

Defining af 5 = (r^-r^) andb| fc) = (rf-r^), Eq. (T£3) expresses that £ fc m fc ((af } ) 2 -(bf } ) 2 ) = £ fe m fc (af } + 
bf')(a| fc) - bf } ) = 0, and thus via Eq. (pS| 

EM" Wrh k x4 k) ) (rffi - r£\) = 0. (A5) 

fc 

Equivalently, 

E^MS-EWi- (A6) 

fc & 
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Combining Eq. (|A6[) with the discrete NVU algorithm (Eq. (|A4p ) gives the following result 



2E fc ^-(r! fc) -r^) 



(A7) 




The atomic NVU algorithm with varying masses is thus given by 




(A8) 



(A9) 
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